---
title: "conjoint analysis Palestine"
author: "Dai Yamao"
date: "`r format(Sys.time(), '%Y-%m-%d')`"
output:
  html_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.width = 10, fig.height = 5)
rm(list=ls())

require(tidyverse)
require(cjoint)
library(cregg)
library(scales)
library(patchwork)
library(extrafont)
require(memisc)
require(stringi)
require(ggplot2)
require(ggpubr)


dat <- read_csv("data/cj_data_Palestine.csv")
dat <- na.exclude(dat)

str(dat)
dat[,c(39:47)] <- data.frame(lapply(dat[,c(39:47)],as.factor))
dat[,c(49:51)] <- data.frame(lapply(dat[,c(49:51)],as.factor))

# filer finished sample
#dat <- dat |> 
#  filter(Finished == "1")

# Ceasefire: unexpected event during survey
dat$RecordedDate <- as.Date(stri_match_first_regex(dat$RecordedDate, "\\d{4}-\\d{2}-\\d{2}"))

dat$ceasefire <- NA
dat$ceasefire[dat$RecordedDate <= as.Date("2025-01-16")] <- "Pre"
dat$ceasefire[as.Date("2025-01-17") <= dat$RecordedDate] <- "Post"
dat$ceasefire <- factor(dat$ceasefire)

# task
dat$task <- as.factor(dat$task)

# Gender
dat <- dat |> 
  mutate(Gender = case_when(Gender == 1 ~ "Male",
                            Gender == 2 ~ "Female"))
dat$Gender <- as.factor(dat$Gender)

# Age group
dat <- dat |> 
  mutate(Age_group = case_when(D2 == 1 ~ "18-20",
                               D2 == 2 ~ "21-30",
                               D2 == 3 ~ "31-40",
                               D2 == 4 ~ "41-50",
                               D2 == 5 ~ "51-60"))
dat$Age_group <- as.factor(dat$Age_group)

# Education
dat <- dat |>
  mutate(Education = D3)
dat$Education <- as.factor(dat$Education)

# Income
dat <- dat |>
  mutate(Income = D5)
dat$Income <- as.factor(dat$Income)

# Sect
dat <- dat |> 
  mutate(Sect = case_when(D6 == 1 ~ "Muslim",
                          D6 == 2 ~ "Christian",
                          D6 == 3 ~ "Jewish",
                          D6 == 4 ~ "Others"))
dat$Sect <- as.factor(dat$Sect)

# Islamism
dat <- dat |> 
  mutate(Islamism = case_when(Q1_1 >= 67 ~ "High",
                              Q1_1 <= 66 & Q1_1 >= 34 ~ "Middle",
                              Q1_1 <= 33 ~ "Low"))
dat$Islamism <- as.factor(dat$Islamism)

# Secularism
dat <- dat |> 
  mutate(Secularism = case_when(Q1_2 >= 67 ~ "High",
                                Q1_2 <= 66 & Q1_2 >= 34 ~ "Middle",
                                Q1_2 <= 33 ~ "Low"))
dat$Secularism <- as.factor(dat$Secularism)

# Localism
dat <- dat |> 
  mutate(Localism = case_when(Q1_3 >= 67 ~ "High",
                               Q1_3 <= 66 & Q1_3 >= 34 ~ "Middle",
                               Q1_3 <= 33 ~ "Low"))
dat$Localism <- as.factor(dat$Localism)

# Arabism
dat <- dat |> 
  mutate(Arabism = case_when(Q1_4 >= 67 ~ "High",
                             Q1_4 <= 66 & Q1_4 >= 34 ~ "Middle",
                             Q1_4 <= 33 ~ "Low"))
dat$Arabism <- as.factor(dat$Arabism)

# US
dat <- dat |> 
  mutate(US = case_when(Q2_1 >= 4 ~ "Like",
                        Q2_1 == 3 ~ "Neither",
                        Q2_1 <= 2 ~ "Dislike"))
dat$US <- as.factor(dat$US)

# Iran
dat <- dat |> 
  mutate(Iran = case_when(Q2_2 >= 4 ~ "Like",
                          Q2_2 == 3 ~ "Neither",
                          Q2_2 <= 2 ~ "Dislike"))
dat$Iran <- as.factor(dat$Iran)

# Saudi
dat <- dat |> 
  mutate(Saudi = case_when(Q2_3 >= 4 ~ "Like",
                           Q2_3 == 3 ~ "Neither",
                           Q2_3 <= 2 ~ "Dislike"))
dat$Saudi <- as.factor(dat$Saudi)

# Turkey
dat <- dat |> 
  mutate(Turkey = case_when(Q2_4 >= 4 ~ "Like",
                            Q2_4 == 3 ~ "Neither",
                            Q2_4 <= 2 ~ "Dislike"))
dat$Turkey <- as.factor(dat$Turkey)

# Qatar
dat <- dat |> 
  mutate(Qatar = case_when(Q2_5 >= 4 ~ "Like",
                          Q2_5 == 3 ~ "Neither",
                          Q2_5 <= 2 ~ "Dislike"))
dat$Qatar <- as.factor(dat$Qatar)

# Israel
dat <- dat |> 
  mutate(Israel = case_when(Q2_6 >= 4 ~ "Like",
                            Q2_6 == 3 ~ "Neither",
                            Q2_6 <= 2 ~ "Dislike"))
dat$Israel <- as.factor(dat$Israel)

# Axis
dat <- dat |> 
  mutate(Axis = case_when(Q2_7 >= 4 ~ "Like",
                          Q2_7 == 3 ~ "Neither",
                          Q2_7 <= 2 ~ "Dislike"))
dat$Axis <- as.factor(dat$Axis)

# UN
dat <- dat |> 
  mutate(UN = case_when(Q2_8 >= 4 ~ "Like",
                        Q2_8 == 3 ~ "Neither",
                        Q2_8 <= 2 ~ "Dislike"))
dat$UN <- as.factor(dat$UN)


# Hamas_Jihad supporter
dat <- dat |> 
  mutate(Hamas_Jihad_supporter = (Q4_2 + Q4_3)/2) |>
  mutate(Hamas_Jihad_supporter = ifelse(Hamas_Jihad_supporter >= 7, "Hamas_Jihad Supporter", "Others"))
dat$Hamas_Jihad_supporter <- as.factor(dat$Hamas_Jihad_supporter)

# Patronage
dat <- dat |> 
  mutate(Patronage = case_when(Q4 >= 3 ~ "High",
                               Q4 <= 2 ~ "Low"))
dat$Patronage <- as.factor(dat$Patronage)

# High patronage and Hamas_Jihad_supporter: does not exist
dat <- dat |>
  mutate(Patron_Hamas = ifelse(Hamas_Jihad_supporter == "Hamas_Jihad Supporter" & Patronage == "High", "Patron_Hamas", "Other"))
dat$Patron_Hamas <- as.factor(dat$Patron_Hamas)

# Rely job on state apparatuses
dat <- dat |> 
  mutate(Job_rely_state = ifelse(Q5_1 <= 5, "Rely job on state", "Others"))
dat$Job_rely_state <- as.factor(dat$Job_rely_state)

# Rely economic difficulty on state apparatuses
dat <- dat |> 
  mutate(economy_rely_state = ifelse(Q5_1 <= 5, "Rely economy on state", "Others"))
dat$economy_rely_state <- as.factor(dat$economy_rely_state)

# Rely security on state apparatuses
dat <- dat |> 
  mutate(Security_rely_state = ifelse(Q5_3 <= 5, "Rely security on state", "Others"))
dat$Security_rely_state <- as.factor(dat$Security_rely_state)

# Piety
dat <- dat |> 
  mutate(Piety = case_when(Q6 <= 4 ~ "Low",
                           Q6 >= 5 & Q6 <= 6 ~ "Middle",
                           Q6 >=7 ~ "High"))
dat$Piety <- as.factor(dat$Piety)

# Subjective conflict intensity
dat <- dat |> 
  mutate(Conflict_assessment = case_when(Q14 >= 4 ~ "Worse",
                                         Q14 == 3 ~ "Same",
                                         Q14 <= 2 ~ "Better"))
dat$Conflict_assessment <- as.factor(dat$Conflict_assessment)

# Subjective security perspective
dat <- dat |> 
  mutate(Future_expectation = case_when(Q15 >= 6 ~ "Worse",
                                        Q15 == 5 ~ "Same",
                                        Q15 <= 2 ~ "Better"))
dat$Future_expectation <- as.factor(dat$Future_expectation)

```


# whole model
```{r, fig.height = 5, dpi = 300}
f1 <- selected ~ Activity + Armament + Patron + Service
str(dat)

# MM
whole_mm <- cj(dat, f1, 
               id = id,
               estimate = "mm", 
               h0 = 0.5,
               alpha = 0.1)
whole_mm
plot(whole_mm, vline = 0.5, 
     size = 2)
write_html(whole_mm, file = "result/whole_mm_Palestine.html")

# amce
whole_amce <- cj(data = dat, f1,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 alpha = 0.1)
whole_amce
plot(whole_amce, vline = 0,
     size = 2)

whole_Palestine <- plot(whole_amce, vline = 0, size = 2) +
  ggplot2::labs(title = "Whole result of Palestine")
ggsave(filename = "result/whole_Palestine.png",
       plot = whole_Palestine,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)

write_html(whole_amce, file = "result/whole_amce_Palestine.html")

# plot estimation difference
whole_mm$Estimate <- "MM"
whole_amce$Estimate <- "AMCE"
plot(rbind(whole_mm, whole_amce), feature_headers = FALSE, size = 2)+
  ggplot2::facet_wrap(~ Estimate, ncol = 2L)
```

# Ceasefire
```{r, fig.height = 5, dpi = 300}
# amce
ceasefire_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ ceasefire,
                 alpha = 0.1)
ceasefire_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ ceasefire,
                      alpha = 0.1)
plot(rbind(ceasefire_amce, ceasefire_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
dat$ceasefire <- relevel(dat$ceasefire, "Pre")
ceasefire_mm <- cj(dat,
               selected ~ Activity + Armament + Patron + Service,
               id = id,
               estimate = "mm",
               by = ~ ceasefire,
               alpha = 0.1)
ceasefire_mm_diff <- cj(data = dat,
                    selected ~ Activity + Armament + Patron + Service,
                    id = id,
                    estimate = "mm_diff",
                    by = ~ ceasefire,
                    alpha = 0.1)
plot(rbind(ceasefire_mm, ceasefire_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(ceasefire_mm, group = "ceasefire", size = 2)

# plot
p1 <- plot(ceasefire_mm, group = "ceasefire", size = 2) +
  ggplot2::labs(title = "Palestine")
p2 <- plot(ceasefire_mm_diff, size = 2)+ 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
p3 <- ggarrange(p1,
                p2 + theme(axis.title.y = element_blank()),
                nrow = 1, ncol = 2, align = "hv",
                common.legend = FALSE,
                legend = "bottom")

ggsave(filename = "result/ceasfire_pl.png",
       plot = p3,
       width = 10,
       height = 4,
       units = "in", 
       dpi = 600)

# OLS
write_html(ceasefire_mm, file = "result/ceasefire_mm_pl.html")
write_html(ceasefire_mm_diff, file = "result/ceasefire_mm_diff_pl.html")

# anova
cj_anova(data = dat, 
         selected ~ Activity + Armament + Patron + Service, by = ~ ceasefire)
```

# task
```{r, fig.height = 7, dpi = 300}
# amce
task_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ task)
task_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ task)
plot(task_amce, size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)
plot(rbind(task_amce, task_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# gender
```{r, fig.height = 5, dpi = 300}
# amce
gender_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Gender)
gender_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Gender)
plot(rbind(gender_amce, gender_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mmm
gender_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Gender)
gender_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Gender)
plot(rbind(gender_mm, gender_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# age group: Design has only one primary sampling unit
```{r, fig.height = 7, dpi = 300}
# amce
#Age_group_amce <- cj(data = dat,
#                      selected ~ Activity + Armament + Patron + Service,
#                      id = ~ respondentIndex,
#                      estimate = "amce",
#                      by = ~ Age_group)
#Age_group_amce_diff <- cj(data = dat,
#                      selected ~ Activity + Armament + Patron + Service,
#                      id = ~ respondentIndex,
#                      estimate = "amce_diff",
#                      by = ~ Age_group)
#plot(Age_group_amce, size = 2, vline = 0) +
#  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
#Age_group_mm <- cj(data = dat,
#                 selected ~ Activity + Armament + Patron + Service,
#                 id = ~ respondentIndex,
#                 estimate = "mm",
#                 by = ~ Age_group)
#Age_group_mm_diff <- cj(data = dat,
#                      selected ~ Activity + Armament + Patron + Service,
#                      id = ~ respondentIndex,
#                     estimate = "mm_diff",
#                      by = ~ Age_group)
#plot(Age_group_mm, size = 2, vline = 0) +
#  ggplot2::facet_wrap(~ BY, ncol = 3L)

#plot(Age_group_mm, group = "Age_group", size = 2)
```

# Education: Design has only one primary sampling unit
```{r, fig.height = 7, dpi = 300}
# amce
#Education_amce <- cj(data = dat,
#                 selected ~ Activity + Armament + Patron + Service,
#                 id = ~ respondentIndex,
#                 estimate = "amce",
#                 by = ~ Education)
#Education_amce_diff <- cj(data = dat,
#                      selected ~ Activity + Armament + Patron + Service,
#                      id = ~ respondentIndex,
#                      estimate = "amce_diff",
#                      by = ~ Education)
#plot(Education_amce, size = 2, vline = 0) +
#  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
#Education_mm <- cj(data = dat,
#                 selected ~ Activity + Armament + Patron + Service,
#                 id = ~ respondentIndex,
#                 estimate = "mm",
#                by = ~ Education)
#Education_mm_diff <- cj(data = dat,
#                      selected ~ Activity + Armament + Patron + Service,
#                      id = ~ respondentIndex,
#                      estimate = "mm_diff",
#                      by = ~ Education)
#plot(Education_mm, size = 2, vline = 0) +
#  ggplot2::facet_wrap(~ BY, ncol = 3L)

#plot(Education_mm, group = "Education", size = 2)
```

# Income
```{r, fig.height = 7, dpi = 300}
# amce
Income_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Income)
Income_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Income)
plot(Income_amce, size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Income_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Income)
Income_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Income)
plot(Income_mm, size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Sect: Design has only one primary sampling unit
```{r, fig.height = 7, dpi = 300}
table(dat$Sect)
# amce
#Sect_amce <- cj(data = dat,
#                 selected ~ Activity + Armament + Patron + Service,
#                 id = ~ respondentIndex,
#                 estimate = "amce",
#                 by = ~ Sect)
#Sect_amce_diff <- cj(data = dat,
#                      selected ~ Activity + Armament + Patron + Service,
#                      id = ~ respondentIndex,
#                      estimate = "amce_diff",
#                      by = ~ Sect)
#plot(rbind(Sect_amce, Sect_amce_diff), size = 2, vline = 0) +
#  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
#Sect_mm <- cj(data = dat,
#                selected ~ Activity + Armament + Patron + Service,
#                 id = ~ respondentIndex,
#                 estimate = "mm",
#                 by = ~ Sect)
#Sect_mm_diff <- cj(data = dat,
#                      selected ~ Activity + Armament + Patron + Service,
#                      id = ~ respondentIndex,
#                      estimate = "mm_diff",
#                      by = ~ Sect)
#plot(rbind(Sect_mm, Sect_mm_diff), size = 2, vline = 0) +
#  ggplot2::facet_wrap(~ BY, ncol = 3L)

#plot(Sect_mm, group = "Sect", size = 2)
```

# Islamism
```{r, fig.height = 7, dpi = 300}
# amce
Islamism_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Islamism)
Islamism_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Islamism)
plot(rbind(Islamism_amce, Islamism_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Islamism_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Islamism)
Islamism_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Islamism)
plot(rbind(Islamism_mm, Islamism_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Islamism_mm, group = "Islamism", size = 2)

```

# Secularism
```{r, fig.height = 7, dpi = 300}
# amce
Secularism_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Secularism)
Secularism_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Secularism)
plot(rbind(Secularism_amce, Secularism_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Secularism_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Secularism)
Secularism_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Secularism)
plot(rbind(Secularism_mm, Secularism_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Secularism_mm, group = "Secularism", size = 2)
```

# Localism
```{r, fig.height = 7, dpi = 300}
# amce
Localism_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Localism)
Localism_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Localism)
plot(rbind(Localism_amce, Localism_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Localism_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Localism)
Localism_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Localism)
plot(rbind(Localism_mm, Localism_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Arabism
```{r, fig.height = 7, dpi = 300}
# amce
Arabism_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Arabism)
Arabism_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Arabism)
plot(rbind(Arabism_amce, Arabism_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Arabism_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Arabism)
Arabism_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Arabism)
plot(rbind(Arabism_mm, Arabism_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Veiws toward US
```{r, fig.height = 7, dpi = 300}
table(dat$US)
# amce
US_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ US)
US_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ US)
plot(rbind(US_amce, US_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
US_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ US)
US_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ US)
plot(rbind(US_mm, US_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(US_mm, group = "US", size = 2)
```

# Veiws toward Iran
```{r, fig.height = 5, dpi = 300}
table(dat$Iran)
# amce
Iran_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Iran,
                alpha = 0.1)
Iran_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Iran,
                     alpha = 0.1)
plot(rbind(Iran_amce, Iran_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Iran_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Iran,
              alpha = 0.1)
Iran_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Iran,
                   alpha = 0.1)
plot(rbind(Iran_mm, Iran_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm plot
plot(Iran_mm, group = "Iran", size = 2)
Iran_Palestine <- plot(Iran_mm, group = "Iran", size = 2) +
  ggplot2::labs(title = "Palestine")
ggsave(filename = "result/Iran_Palestine.png",
       plot = Iran_Palestine,
       width = 7,
       height = 5,
       units = "in", 
       dpi = 600)

# mm difference plot
plot(Iran_mm_diff, size = 2, vline = 0) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
Iran_diff_PL <- plot(Iran_mm_diff, size = 2, vline = 0) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L) +
  ggplot2::labs(title = "Palestine")
ggsave(filename = "result/Iran_diff_PL.png",
       plot = Iran_diff_PL,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)

# anova
cj_anova(data = dat, 
         selected ~ Activity + Armament + Patron + Service, by = ~ Iran)

# OLS result
write_html(Iran_mm, file = "result/Iran_mm_pl.html")
write_html(Iran_mm_diff, file = "result/Iran_mm_diff_pl.html")
```

# Veiws toward Saudi
```{r, fig.height = 7, dpi = 300}
table(dat$Saudi)
# amce
Saudi_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Saudi)
Saudi_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Saudi)
plot(rbind(Saudi_amce, Saudi_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Saudi_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Saudi)
Saudi_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Saudi)
plot(rbind(Saudi_mm, Saudi_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Veiws toward Turkey
```{r, fig.height = 7, dpi = 300}
table(dat$Turkey)
# amce
Turkey_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Turkey)
Turkey_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Turkey)
plot(rbind(Turkey_amce, Turkey_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Turkey_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Turkey)
Turkey_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Turkey)
plot(rbind(Turkey_mm, Turkey_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Veiws toward Qatar
```{r, fig.height = 7, dpi = 300}
table(dat$Qatar)
# amce
Qatar_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Qatar)
Qatar_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Qatar)
plot(rbind(Qatar_amce, Qatar_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Qatar_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Qatar)
Qatar_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Qatar)
plot(rbind(Qatar_mm, Qatar_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Veiws toward Israel
```{r, fig.height = 7, dpi = 300}
table(dat$Israel)
# amce
Israel_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Israel)
Israel_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Israel)
plot(rbind(Israel_amce, Israel_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Israel_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Israel)
Israel_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Israel)
plot(rbind(Israel_mm, Israel_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Veiws toward Axis of resistance
```{r, fig.height = 7, dpi = 300}
table(dat$Axis)
# amce
Axis_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Axis)
Axis_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Axis)
plot(rbind(Axis_amce, Axis_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Axis_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Axis)
Axis_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Axis)
plot(rbind(Axis_mm, Axis_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Axis_mm, group = "Axis", size = 2)
```

# Veiws toward UN
```{r, fig.height = 7, dpi = 300}
table(dat$UN)
# amce
UN_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ UN)
UN_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ UN)
plot(rbind(UN_amce, UN_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
UN_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ UN)
UN_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ UN)
plot(rbind(UN_mm, UN_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Hamas and Islamic Jihad supporter
```{r, fig.height = 5, dpi = 300}
table(dat$Hamas_Jihad_supporter)
# amce
Hamas_Jihad_supporter_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Hamas_Jihad_supporter,
                 alpha = 0.1)
Hamas_Jihad_supporter_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Hamas_Jihad_supporter,
                      alpha = 0.1)
plot(rbind(Hamas_Jihad_supporter_amce, Hamas_Jihad_supporter_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Hamas_Jihad_supporter_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Hamas_Jihad_supporter,
                 alpha = 0.1)
Hamas_Jihad_supporter_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Hamas_Jihad_supporter,
                      alpha = 0.1)
plot(rbind(Hamas_Jihad_supporter_mm, Hamas_Jihad_supporter_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Hamas_Jihad_supporter_mm, group = "Hamas_Jihad_supporter", size = 2)

```

# Patronage
```{r, fig.height = 5, dpi = 300}
# amce
Patronage_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Patronage,
                 alpha = 0.1)
Patronage_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Patronage,
                      alpha = 0.1)
plot(rbind(Patronage_amce, Patronage_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Patronage_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Patronage,
                 alpha = 0.1)
Patronage_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Patronage,
                      alpha = 0.1)
plot(rbind(Patronage_mm, Patronage_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Patronage_mm, group = "Patronage", size = 2)
```

# High patronage and Hamas Jihad supporter: Design has only one primary sampling unit
```{r, fig.height = 5, dpi = 300}
# amce
#Patron_Hamas_amce <- cj(data = dat,
#                 selected ~ Activity + Armament + Patron + Service,
#                 id = ~ respondentIndex,
#                 estimate = "amce",
#                 by = ~ Patron_Hamas)
#Patron_Hamas_amce_diff <- cj(data = dat,
#                      selected ~ Activity + Armament + Patron + Service,
#                      id = ~ respondentIndex,
#                      estimate = "amce_diff",
#                      by = ~ Patron_Hamas)
#plot(rbind(Patron_Hamas_amce, Patron_Hamas_amce_diff), size = 2) + 
#  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Rely job on state apparatuses
```{r, fig.height = 5, dpi = 300}
table(dat$Job_rely_state)
# amce
Job_rely_state_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Job_rely_state)
Job_rely_state_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Job_rely_state)
plot(rbind(Job_rely_state_amce, Job_rely_state_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Job_rely_state_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Job_rely_state)
Job_rely_state_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Job_rely_state)
plot(rbind(Job_rely_state_mm, Job_rely_state_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```


# Rely economic dificulty on state apparatuses
```{r, fig.height = 5, dpi = 300}
table(dat$economy_rely_state)
# amce
economy_rely_state_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ economy_rely_state)
economy_rely_state_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ economy_rely_state)
plot(rbind(economy_rely_state_amce, economy_rely_state_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
economy_rely_state_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ economy_rely_state)
economy_rely_state_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ economy_rely_state)
plot(rbind(economy_rely_state_mm, economy_rely_state_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Rely security on state apparatuses
```{r, fig.height = 5, dpi = 300}
table(dat$Security_rely_state)
# amce
Security_rely_state_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Security_rely_state)
Security_rely_state_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Security_rely_state)
plot(rbind(Security_rely_state_amce, Security_rely_state_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Security_rely_state_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Security_rely_state)
Security_rely_state_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Security_rely_state)
plot(rbind(Security_rely_state_mm, Security_rely_state_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Piety
```{r, fig.height = 7, dpi = 300}
table(dat$Piety)
# amce
Piety_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Piety)
Piety_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Piety)
plot(rbind(Piety_amce, Piety_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Piety_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Piety)
Piety_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Piety)
plot(rbind(Piety_mm, Piety_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Piety_mm, group = "Piety", size = 2)
```

# Conflict_assessment: Design has only one primary sampling unit
```{r, fig.height = 7, dpi = 300}
table(dat$Conflict_assessment)
# amce
#Conflict_assessment_amce <- cj(data = dat,
#                 selected ~ Activity + Armament + Patron + Service,
#                 id = ~ respondentIndex,
#                 estimate = "amce",
#                 by = ~ Conflict_assessment)
#Conflict_assessment_amce_diff <- cj(data = dat,
#                      selected ~ Activity + Armament + Patron + Service,
#                      id = ~ respondentIndex,
#                      estimate = "amce_diff",
#                      by = ~ Conflict_assessment)
#plot(rbind(Conflict_assessment_amce, Conflict_assessment_amce_diff), size = 2, vline = 0) +
#  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
#Conflict_assessment_mm <- cj(data = dat,
#                 selected ~ Activity + Armament + Patron + Service,
#                 id = ~ respondentIndex,
#                 estimate = "mm",
#                 by = ~ Conflict_assessment)
#Conflict_assessment_mm_diff <- cj(data = dat,
#                      selected ~ Activity + Armament + Patron + Service,
#                      id = ~ respondentIndex,
#                      estimate = "mm_diff",
#                      by = ~ Conflict_assessment)
#plot(rbind(Conflict_assessment_mm, Conflict_assessment_mm_diff), size = 2, vline = 0) +
#  ggplot2::facet_wrap(~ BY, ncol = 3L)

#plot(Conflict_assessment_mm, group = "Conflict_assessment", size = 2)
```

# Future_expectation
```{r, fig.height = 7, dpi = 300}
table(dat$Future_expectation)
# amce
Future_expectation_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Future_expectation)
Future_expectation_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Future_expectation)
plot(rbind(Future_expectation_amce, Future_expectation_amce_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Future_expectation_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Future_expectation)
Future_expectation_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Future_expectation)
plot(rbind(Future_expectation_mm, Future_expectation_mm_diff), size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Future_expectation_mm, group = "Future_expectation", size = 2)
```

# respondant location
```{r}
table(dat$LocationLatitude, dat$LocationLongitude)
```

# frequency and propotion
```{r}
plot(cj_freqs(dat,
              f1,
              id = ~ respondentIndex))
```

# interaction: Ceasefire*Hamas supporter
```{r, fig.height = 5, dpi = 300}
dat$Hamas_ceasefire <- interaction(dat$Hamas_Jihad_supporter, dat$ceasefire, sep = "_")
Hamas_ceasefire_mm <- cj(dat,
               selected ~ Hamas_ceasefire,
               id = id,
               estimate = "mm",
               alpha = 0.1)
plot(Hamas_ceasefire_mm, vline = 0.5, size = 2)

# plot
Hamas_ceasefire <- plot(Hamas_ceasefire_mm, vline = 0, size = 2) +
  labs(title = "Palestine")
ggsave(filename = "result/Hamas_ceasefire.png",
       plot = Hamas_ceasefire,
       width = 5,
       height = 3,
       units = "in", 
       dpi = 600)

# OLS result
write_html(Hamas_ceasefire_mm, file = "result/Hamas_ceasefire_m.html")
```

# interaction
```{r}
table(dat$Iran, dat$Hamas_Jihad_supporter)
cor(as.numeric(dat$Iran), as.numeric(dat$Hamas_Jihad_supporter))

# supporter or not
dat$ceasefire <- relevel(dat$ceasefire, "Pre")
supporter <- dat |>
  filter(Hamas_Jihad_supporter == "Hamas_Jihad Supporter")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm", 
     by = ~ ceasefire,
     alpha = 0.1)
plot(supporter, group = "ceasefire", 
     size = 2,
     vline = 0.5)

non_supporter <- dat |>
  filter(Hamas_Jihad_supporter == "Others")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm", 
     by = ~ ceasefire,
     alpha = 0.1)
plot(non_supporter, group = "ceasefire", 
     size = 2,
     vline = 0.5)+
  aes(shape = ceasefire)+
  scale_shape_manual(values = c(1,2,3,4,5))+
  scale_colour_manual(values=rep("black", 5)) 


plot(rbind(supporter, non_supporter), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# difference
supporter_diff <- dat |>
  filter(Hamas_Jihad_supporter == "Hamas_Jihad Supporter")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm_diff", 
     by = ~ ceasefire,
     alpha = 0.1)
plot(supporter_diff, group = "ceasefire", 
     size = 2)

non_supporter_diff <- dat |>
  filter(Hamas_Jihad_supporter == "Others")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm_diff", 
     by = ~ ceasefire,
     alpha = 0.1)
plot(non_supporter_diff, group = "ceasefire", 
     size = 2)


plot(rbind(supporter_diff, non_supporter_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)


# interaction
supporter$Supporter <- "Supporter"
non_supporter$Supporter <- "Non_supporter"
interaction <- rbind(supporter, non_supporter)

H2_app_pl <- plot(interaction, group = "Supporter", size = 2) +
  ggplot2::facet_wrap(~ BY, ncol = 2) +
  labs(title = "Palestine",
       shape = "Hamas Supporter",
       labels = c("Non-suporter", "Supporter"))+
  scale_color_hue(name = "Hamas Supporter", labels = c("Non-suporter", "Supporter"))
H2_app_pl

ggsave(filename = "result/H2_app_pl.png",
       plot = H2_app_pl,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)  

# OLS result
write_html(interaction, file = "result/H2_app_pl.html")

# diff
supporter_diff$Supporter <- "Supporter"
non_supporter_diff$Supporter <- "Non_supporter"
interaction_diff <- rbind(supporter_diff, non_supporter_diff)

H1_app_pl <- plot(interaction_diff, group = "Supporter", size = 2) +
  ggplot2::facet_wrap(~ BY, ncol = 2) +
  labs(title = "Palestine",
       labels = c("Non-suporter", "Supporter"))+
  scale_color_hue(labels = c("Non-suporter", "Supporter"))
H1_app_pl

ggsave(filename = "result/H1_app_pl.png",
       plot = H1_app_pl,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)  
```
